# Clean up
rm(list = ls())

# Set working directory; please set your own here
setwd("~/Dropbox/Fragmentation/fragmentation_replication_bjps/")

library(tidyverse)
library(readstata13)
library(gridExtra)
library(foreign)

# Import the election-level dataset
first_stage <- read.dta13("Data/fragmentation_electionlevel.dta", generate.factors = TRUE)
data <- read.dta("Results/analyses_final.dta")

# Create a theme for the plots
theme_base <- 
  theme_minimal(base_size = 9)  + 
  theme(legend.position = c(0.85, 0.1),
        axis.text = element_text(size = 12),
        axis.title.x=element_text(size = 12),
        axis.title.y=element_text(size = 12),
        plot.title = element_text(size = 12, hjust = 0.5),
        legend.text = element_text(size = 12))

# Plotting ENPP against ENPP added by parties in the bandwidth
enpp_plot <- ggplot(first_stage, aes(x = enppadd_bw50, y = enpp)) +
  geom_hline(yintercept = 0, linetype = "dashed", color = "red",size = .5) +
  geom_point(alpha = 0.6) +
  stat_smooth(method = "lm", aes(colour = "Linear model")) +
  stat_smooth(method = "loess", linetype = "longdash", aes(colour = "LOESS")) +
  scale_colour_manual(name = " ", values = c("blue", "green3")) +
  scale_y_continuous(limits = c(-3, 10)) +
  theme_base +
  labs(x = "Effective number of parl parties added by parties within the bandwidth", y = "Effective number of parliamentary parties")

# Plot distributuion of f-statistics in first stage
firststage_distr_enpp <- data %>%
  filter(optimal == 1) %>%
  filter(Sample == "Whole") %>%
  filter(FE == "None") %>%
  filter(id == "No. parties just above") %>%
  ggplot +
  scale_x_continuous(name = "F-statistic", limits = c(20, 37.5)) +
  scale_y_continuous(name = "Density") +
  geom_vline(aes(xintercept = mean(fstat)), col = 'red', size = 1, linetype = "dashed") +
  theme_bw() +
  geom_density(aes (x = fstat), fill = "steelblue", color = "blue", alpha = 0.5) 
firststage_distr_enpp

# Check mean f-stat in the main specification
data_2 <- data %>%
  filter(optimal == 1) %>%
  filter(Sample == "Whole") %>%
  filter(FE == "None") %>%
  filter(id == "No. parties just above")

mean(data_2$fstat)

firststage_main <- grid.arrange(enpp_plot, firststage_distr_enpp, ncol = 2, nrow = 1)
firststage_main

# Save the final plot
ggsave("Plots/fig3.png", plot = firststage_main, width = 30, height = 20, units = "cm")